#include <iostream>
#include <fstream>
#include <sstream>
#include "ParmParse.H"
#include "LevelData.H"
#include "FArrayBox.H"
#include "Vector.H"
#include "aceMesh_runtime.h"

static Box
boxFromVector(const Vector<int>&  a_ints,
              int                 a_start=0)
{
  IntVect lo(a_ints[a_start+0], a_ints[a_start+1], a_ints[a_start+2]);
  IntVect hi(a_ints[a_start+3], a_ints[a_start+4], a_ints[a_start+5]);
  Box bx(lo, hi);
  return bx;
}

static void 
initialize_grids(DisjointBoxLayout &a_dbl)
{
  //Read boxes from inputfiles
  ParmParse pp;
  vector<int> nCells(SpaceDim);
  pp.getarr("n_cells",nCells,0,SpaceDim);
  Real coarseDx = 1.0/nCells[0];

  IntVect lo = IntVect::Zero;
  IntVect hi;
  for (int idir = 0 ; idir < SpaceDim ; idir++){
    hi[idir] = nCells[idir];
  }
  Box DomBox(lo,hi);
  ProblemDomain proDom(DomBox);

  Vector<Box> boxes;
  int level1boxcount;
  Vector<int> level1domainIndices,level1boxIndices;
  Vector<int> procAssign;

  if (pp.query("level_1_boxcount", level1boxcount))
  {
    pp.getarr("level_1_boxes", level1boxIndices, 0, 2*3*level1boxcount);
  }

  for (int indbox = 0; indbox < level1boxcount; indbox++)
  {
    Box baseBox = boxFromVector(level1boxIndices, 2*3*indbox);
    boxes.push_back(baseBox);
    procAssign.push_back(0);
  }
  a_dbl.define(boxes, procAssign, proDom);
}


static void 
run_cases(LevelData<FArrayBox> &a_lhs,
          LevelData<FArrayBox> &a_rhs)
{
  int numcomp = a_lhs.nComp();
  int startcomp = 0;
  DataIterator dit = a_lhs.dataIterator();
//#pragma acemesh for 
  for(dit.begin(); dit.ok(); ++dit)
  { 
    DataIndex d = dit();
    Box subbox = a_lhs.disjointBoxLayout()[d];
    //a_lhs[d].plus(a_rhs[d], subbox, subbox, 1.0, startcomp, startcomp, numcomp);
	
    FArrayBox& destfab = a_lhs[d];
    FArrayBox& srcfab = a_rhs[d];

//#pragma acemesh data(destfab:W,srcfab:R)	
    if (!((subbox).isEmpty())) 
    {                                             
      const int *_th_plo = (destfab).loVect();                            
      IntVect _th_plen = (destfab).size();                             
      const int *_x_plo = (srcfab).loVect();                               
      IntVect _x_plen = (srcfab).size();                                
      const int *_subbox_lo = (subbox).loVect();                           
      IntVect _subbox_len = (subbox).size();                  
      const int *_bx_lo = (subbox).loVect();                              
      Real* _th_p = (destfab).dataPtr(0);                                   
      const Real* _x_p  = (srcfab).dataPtr(0);                             

      for (int _k = 0; _k < _subbox_len[2]; ++_k) {               
        const int kR = _k + _subbox_lo[2];                      
        const int kxR = _k + _bx_lo[2];                     
        for (int _j = 0; _j < _subbox_len[1]; ++_j) {           
          const int jR = _j + _subbox_lo[1];                  
          const int jxR = _j + _bx_lo[1];                 
          Real *_th_pp = _th_p + ((_subbox_lo[0] - _th_plo[0])                 
              + _th_plen[0]*((jR - _th_plo[1])
                + _th_plen[1]*(kR - _th_plo[2])));               
          const Real *_x_pp = _x_p + ((_bx_lo[0] - _x_plo[0])                      
              + _x_plen[0]*((jxR - _x_plo[1])                    
                + _x_plen[1]*(kxR - _x_plo[2])));                
          for (int _i = 0; _i < _subbox_len[0]; ++_i, ++_th_pp) {                              
            int iR = _i + _subbox_lo[0]; 
            Real &thisR = * _th_pp; 
            const Real & xR = _x_pp[_i];
            thisR += xR * 1.0;
          }
        }
      }
    }
	
  }
}

typedef struct 
{
	FArrayBox *_lhs;
	FArrayBox *_rhs;
	const Box *_subbox;
}targs;

void _do_job(targs *args)
{
	if (!((*(args->_subbox)).isEmpty())) 
    {                                             
      const int *_th_plo = (*(args->_lhs)).loVect();                            
      IntVect _th_plen = (*(args->_lhs)).size();                             
      const int *_x_plo = (*(args->_rhs)).loVect();                               
      IntVect _x_plen = (*(args->_rhs)).size();                                
      const int *_subbox_lo = (*(args->_subbox)).loVect();                           
      IntVect _subbox_len = (*(args->_subbox)).size();                  
      const int *_bx_lo = (*(args->_subbox)).loVect();                              
      Real* _th_p = (*(args->_lhs)).dataPtr(0);                                   
      const Real* _x_p  = (*(args->_rhs)).dataPtr(0);                             

      for (int _k = 0; _k < _subbox_len[2]; ++_k) {               
        const int kR = _k + _subbox_lo[2];                      
        const int kxR = _k + _bx_lo[2];                     
        for (int _j = 0; _j < _subbox_len[1]; ++_j) {           
          const int jR = _j + _subbox_lo[1];                  
          const int jxR = _j + _bx_lo[1];                 
          Real *_th_pp = _th_p + ((_subbox_lo[0] - _th_plo[0])                 
              + _th_plen[0]*((jR - _th_plo[1])
                + _th_plen[1]*(kR - _th_plo[2])));               
          const Real *_x_pp = _x_p + ((_bx_lo[0] - _x_plo[0])                      
              + _x_plen[0]*((jxR - _x_plo[1])                    
                + _x_plen[1]*(kxR - _x_plo[2])));                
          for (int _i = 0; _i < _subbox_len[0]; ++_i, ++_th_pp) {                              
            int iR = _i + _subbox_lo[0]; 
            Real &thisR = * _th_pp; 
            const Real & xR = _x_pp[_i];
            thisR += xR * 1.0;
    }}}}
}


static void 
run_cases_cdag(LevelData<FArrayBox> &a_lhs,
          LevelData<FArrayBox> &a_rhs)
{
  int numcomp = a_lhs.nComp();
  int startcomp = 0;
  DataIterator dit = a_lhs.dataIterator();
  
  begin_split_task();
  for(dit.begin(); dit.ok(); ++dit)
  { 
    DataIndex d = dit();
    const Box& subbox = a_lhs.disjointBoxLayout()[d];
    //a_lhs[d].plus(a_rhs[d], subbox, subbox, 1.0, startcomp, startcomp, numcomp);
	
    FArrayBox& destfab = a_lhs[d];
    FArrayBox& srcfab = a_rhs[d];

	targs args;
	args._lhs = &destfab;
	args._rhs = &srcfab;
	args._subbox = &subbox;
	
	__acemesh_push_wrlist(1, &destfab, UNSHADE);
	__acemesh_push_rlist(1, &srcfab, UNSHADE);
	__acemesh_task_generator((TASK_FUNCPTR)_do_job, &args, sizeof(targs));
  }
  end_split_task();
  spawn_and_wait();
}

static void
outputData(LevelData<FArrayBox> &lhs, const char* filename)
{
  ofstream fout(filename, ios::app);
  DataIterator dit = lhs.dataIterator();
  for (dit.begin(); dit.ok(); ++dit)
  {
    DataIndex d = dit();
    FArrayBox& fab = lhs[d];
    Real* _th_p = fab.dataPtr(0);
    const int* _th_plo = fab.loVect();
    IntVect _th_plen = fab.size();
    const Box &b = fab.box();
    const int* _b_lo = b.loVect();
    IntVect _b_len = b.size();

    fout<<d.intCode()<<" "<<b.smallEnd()<<","<<b.bigEnd()<<endl;
    for(int k=0; k<_b_len[2]; k++){
      const int kR = k + _b_lo[2];
      for(int j =0; j<_b_len[1]; j++){
        const int jR = j + _b_lo[1];
        Real *_th_pp = _th_p + ((_b_lo[0]-_th_plo[0])+
            _th_plen[0]*((jR-_th_plo[1])+
              _th_plen[1]*(kR-_th_plo[2])));
        for (int i=0; i<_b_len[0];++i, ++_th_pp){
          fout<<*_th_pp<<" ";
        }
        fout<<endl;
      }
    }
    fout<<endl;
  }

  fout.close();
}


int main(int argc, char* argv[])
{
  if (argc < 2)
  {
    cerr<< " usage " << argv[0] << " <input_file_name> " << endl;
    exit(0);
  }
  char* in_file = argv[1];
  ParmParse  pp(0,argv+2,NULL,in_file);

  DisjointBoxLayout dbl;
  LevelData<FArrayBox> *lhs;
  LevelData<FArrayBox> *lhs_dag;
  LevelData<FArrayBox> *rhs;

  initialize_grids(dbl);

  lhs = new LevelData<FArrayBox>(dbl, 1, 2*IntVect::Unit);
  lhs_dag = new LevelData<FArrayBox>(dbl, 1, 2*IntVect::Unit);
  rhs = new LevelData<FArrayBox>(dbl, 1, 2*IntVect::Unit);
  for(DataIterator dit = dbl.dataIterator(); dit.ok(); ++dit)
  {
    (*lhs)[dit()].setVal(0.);
	(*lhs_dag)[dit()].setVal(0.);
    (*rhs)[dit()].setVal(1.);
  }

  run_cases(*lhs, *rhs);
  run_cases_cdag(*lhs_dag, *rhs);

  if (argc == 3)
  {
    char *out_file = argv[2];
    outputData(*lhs, out_file);
	
	stringstream filename;
	filename<<out_file;
	filename<<"_dag";

	outputData(*lhs_dag, filename.str().c_str());
  }

  pp.undefine();

  return 0;
}


